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ABSTRACT 

Ionization front instabilities have long been of interest for their suspected role in a variety of phe- 
nomena in the galaxy, from the formation of bright rims and 'elephant trunks' in nebulae to triggered 
star formation in molecular clouds. Numerical treatments of these instabilities have historically been 
limited in both dimensionality and input physics, leaving important questions about their true evo- 
lution unanswered. We present the first three-dimensional radiation hydrodynamical calculations of 
both R-type and D-type ionization front instabilities in galactic environments (i.e., solar metallicity 
gas). Consistent with linear stability analyses of planar D-type fronts, our models exhibit many short- 
wavelength perturbations growing at early times that later evolve into fewer large-wavelength struc- 
tures. The simulations demonstrate that both self-consistent radiative transfer and three-dimensional 
flow introduce significant morphological differences to unstable modes when compared to earlier two- 
dimensional approximate models. We find that the amplitude of the instabilities in the nonlinear 
regime is primarily determined by the efficiency of cooling within the shocked neutral shell. Strong 
radiative cooling leads to long, extended structures with pronounced clumping while weaker cooling 
leads to saturated modes that devolve into turbulent flows. These results suggest that expanding H 
II regions may either promote or provide turbulent support against the formation of later generations 
of stars, with potential consequences for star formation rates in the galaxy today. 

Subject headings: H II regions: simulation — interstellar medium: ionization fronts — galaxy 



1. INTRODUCTION 

This is the first of two papers that study ionization 
front (I-front) instabilities in the nonlinear regime us- 
ing three-dimensional radiation hydrodynamical simula- 
tions. In this paper we consider I-fronts propagating in 
a gas of solar metallicity, appropriate to galactic mas- 
sive star forming environments. In a second paper, we 
simulate I-front propagation in primordial gas with low 
metallicities relevant to star formation at high redshifts. 
We demonstrate that the cooling properties of the gas 
have a large infiuence on how I-front instabilities mani- 
fest. 

These instabilities have been the subject of numer- 
ous analytical and numerical studies for the past fifty 
years, primarily for their suspected role in the for- 
mation of inhomogeneities obser ved in galac ti c emis- 
sion nebulae since the 1930's (iKah^ flOfiSL lAxfordl 
1964 iGarcia -Seg ura fc Francol[l996t [Franco et al.lll99a 
WiUiamsl 11999. 2002). Initial corrugations in the fronts 
due to ambient density fiuctuations are thought to elon- 
gate into ionized fingers extending ahead of the front that 
can evolve into the 'elephant trunks', bright rims, and 
cometary globules commonly seen in the galaxy today. 
Dynamical instabilities in H II regions are of direct rele- 
vance to the environments of massive stars (|Frever et al.l 
[2003) and gamma-ray burst progenitors in general, possi- 
bly acting in concert with stellar winds to drive clumping 
of gas that later become susceptible to gravitational in- 
stabilities and collapse. Perturbations in the ionization 
front shocks of OB associations may also come into play 
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in scen arios of triggered star formation within molecular 
clouds (|Dale et al.ll2007t ). Historically, studies of I-front 
instabilities have progressed in both dimensionality and 
physics along two lines: u nstable modes in R -type fronts 
and in D-type I-fronts (see lOsterbrocld ( |1989fl for a review 
of mo dern I-front nomenclature and lWhalen fc Normanl 
(|2006l ) for numerical examples in one dimension). 

1.1. Shadowing Instabilities in R-Type Fronts 

R-type I-fronts propagate supersonically through me- 
dia with minimal hydrodynamical coupling to the 
gas. The stability of this type of front to encoun- 
ter s with modest dens ity p erturbations was ex amined 
bv IVandervoortI (I1962D a nd IN ewman fc AxfordI (ji967). 
iNewman fc AxfordI ()1967f ) found the fronts to be weakly 
unstable but were unable to follo w their evolution into 
D-type structures. iBerto ldi' (1989) studied the flash ion- 
ization of diffuse clumps in which the front does not re- 
vert to D-type (the 'cloud-zapping regime'). He con- 
cluded that the clumps evaporate without forming new 
inhomogeneities in their wake but could not determine 
the downstream fate of the front itself. Several groups 
have in vestigated the shadowing of I-fronts by dense 
clum ps (jMellema et al.l Il998t ICa.nto et all 119981 : ISokeil 
|1998[) , with particular focus on the evolution of the neu- 
tral tails behi nd the clumps. Most recent in this vein are 
the studies by 'Shap iro et al.l (|2004 !) of cosmological mini- 
halo photoevaporation by external Population HI and 
quasar ionizing sources at high and intermediate red- 
shifts. The aim of these simulations was to assess the 
role of minihalos as sinks of ionizing UV photons during 
cosmological reionization. Shadows rather than insta- 
bilities arose in these models because the obstructions 
tr apped the I-fron ts as they evaporated. 

IWilliami ()1999[ ) revisited instability modes in R-type 
fronts, in this instance following their transition to D- 
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type far downstream of the original perturbation. Their 
analysis revealed that the overdensity indents a planar 
R-type front with a dimple that fir st elongates but then 
saturates (Fig 1 in lWilliamd (|1999[ )V The ionizing pho- 
tons are oblique to the sides of the dimple but incident to 
its tip. Since the flux is lower at the sides than at the tip, 
as the front approaches the Stromgren position the shock 
breaks through the front along the wall of the dimple be- 
fore its tip. Still R-type, the tip races forward, elongat- 
ing the dimple into a nonlinear jet instability when the 
transition to D-critical is complete. If the shocked gas 
efficiently cools, the thin-shell hydrodynamic instability 
appears and fragments the gas into clumps with densities 
hundreds of times greater than the ambient medium. 

1.2. Dynamical Instabilities in D-Type Fronts 

D-type I-fronts move supersonically with respect to 
the upstream neutral gas, but subsonically (or sonically) 
with respect to the downstream ionized gas. Conse- 
quently, D-type fronts are preceded by shocks in the 
neutral gas. Soon afte r the moder n no menclature for 
I-fronts was established, iKahnI (|1958f) and|A2cford| (|1964f ) 
examined the stability of perturbations in planar weak 
D-type fronts. They discovered that recombinations to 
the ground state within ionized fingers advancing ahead 
of the front effectively dampened their further growth 
upon reaching lengths of a few ni~^ pc, where n^ is the 
electron density in the postfront gas. This is the product 
of the sound speed of the ionized gas and its recombi- 
nation time; in galactic environments such corrugations 
are initially short wavelength and saturate quickly. As 
the fingers lengthen, photons must cross more neutrals 
resulting from recombinations, which attenuate the ion- 
izing flux advancing the tip of the disturbance. 

Instabilities that do exhibit strong growth were later 
discovered in D-type ionization fronts, but onl y in con- 
junct i on with pre-existing flow i nstabilities (jVishniad 
119831 ICapriottilll973l : lBrandlll981f ). In general, expand- 
ing spherical shocks are prone to a variety of dynamical 
instabilities that can result in their fragmentation and 
breakup. Radiating shocks are subject to the Vishniac 
(or thin-shell) instability in either uniform densities or 
gradients, while Rayleigh- Taylor instabilities can arise in 
adiabatic or radiating flows if they accelerate (which oc- 
curs in radial power-law density gradients steeper than 

The Vi shniac mechanism can be vi sualized as follows 
(Fig 1 of lMac Low fc NormanI ([Till)): in the frame of 
the shock, the inflow ram pressure is unidirectional while 
the ionized postshock pressure is isotropic. If a density 
fluctuation is advected across the shock, small transverse 
velocities arise across the face of the shock that create 
overdensities and underdensities along adjacent lines of 
sight from the center of the shell. The postshock gas then 
protrudes through the underdense regions, dimpling the 
shock into peaks and valleys and forcing even greater 
transverse flows that further bunch gas across adjoin- 
ing lines of sight. However, the higher pressures in the 
troughs eventually reverse the transverse flows, resulting 
in oscillatory ripples with modest growth rates across the 
face of the shock commonly known as overstabilities. 

Both thin-shell overstabilities and Rayleigh- Taylor in- 
stabilities are exacerbat e d when driven by an ioniza- 
tion front (jSpitzeil 119541 : iFriemanI 11954 iGiulianil 119791 : 
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pear in the shell, the radiation preferentially advances 
into the lower optical depths of the underdensities, re- 
inforcing the transverse velocity gradients that origi- 
nally created them in a repeat of the cycle described 
above (Fig 1 of Garcia-Segura fc Franco). What dis- 
tinguishes unstable modes in blast waves or wind-blown 
bubbles from those in D-type fronts is their ampli- 
tude: in extreme cases the radiation may escape through 
cracks in th e shock, driving even greater hydrodynami- 
cal g r owth (iGarcia-Segura fc Francol Il996l : iFranco et al.l 
119981 : iFrever et al.ll2003l ). The violence of flow instabili- 
ties in general is largely governed by the cooling rates of 
the shocked gas. Efflcient cooling collapses the shocked 
neutral gas into a thin, dense cold shell more prone to 
breakup than a thick semi-adiabatic shock. We note 
that a class of long- wavelength unstable surface modes in 
weak D-type fronts has recently bee n uncov e red by an- 
alytical a nd numerical work done by ISvsoevI (|1997f l and 
IWilliami |2002.) that are independent of any other flow 
instabilities, but we do not consider them here. 

A key question is how the dimensionality of ionized 
flows dictates the nature of perturbation growth in R- 
type and D-type fronts, given that published simulations 
have been confined to two dimensions. Furthermore, 
instability growth is mediated by I-front transitions in 
prefront shocked gas that mandate radiation hydrody- 
namics rather than the approximations heretofore per- 
formcd (Garcia-Segura fc Francolll996l : iDale et al.ll2007t 
Yoshida et al. 2006). In order to explore these issues we 
have performed the first three-dimensional radiation hy- 
drodynamical calculatio ns of dynamical instabilities in 
D-type fronts (jFranco et al. 199^ and s hadowing insta- 
bilities in R-type fronts (Williamg |1999f ) with metal-line 
cooling. These models have been made possible by re- 
cent upgrades to the ZEUS -MP reactive flow physics code 
(jWhalen fc NormanI [2001 ) (hereafter WN06) to perform 
fully-parallelized three-dimensional radiation transport 
for a point source in spherical polar grids and for plane 
waves in cartesian boxes. Our calculations are intended 
to survey the qualitative nature of three-dimensional in- 
stabilities in galactic ionization fronts rather than pro- 
vide detailed comparisons to analytical studies, which 
have only been performed in one dimension. We focus on 
dynamical phenomena in r^^ density fields for compar- 
ison to earlier two-dimensional numerical work and for 
their similarity to clump profiles observed in massive star 
form ing regions within molecular clouds in the galaxy to- 
davJArauilla fc Goldsmithl[l985l IC^ gorio Hetem et al.l 
Il988f l. Our numerical suite is therefore of direct rele- 
vance to the environments of massive stars as well as 
the morphologies of emission nebulae in the interstellar 
medium (ISM). Since our simulations all employ radia- 
tive cooling in uniform or r~^ densities, we expect initial 
perturbation growth to be dominated by Vishniac modes. 
In nearly adiabatic primordial gas with steeper gradients 
such as those in cosmological minihalos, unstable modes 
could be initiated through the Rayleigh- Taylor channel, 
a possibility we investigate in a following paper. 

In Section 2 we describe the improvements to our ra- 
diative transfer and reaction network integration scheme 
required by transport in three dimensions as well as 
the parallelization and load balancing of the transport 
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itself. In section 3 we compare two-dimensional cal- 
culations of dynamical instabilities on polar grids per- 
formed with simplified ionization equilibrium physics to 
those with full radiative transfer. Our three-dimensional 
simulations of these instabilities on spherical-polar co- 
ordinate grids are detailed in Section 4 and we present 
three-dimension al calculations of the shadow instability 
IWilliami (119991 ) in section 5. Implications of this work 
are discussed in Section 6. 

2. ZEUS-MP ALGORITHM UPGRADES 

The ZEUS-MP hydrocode utihzed in previous one- 
dimensional work (jWhalen et al. Il2004[ ) is now fully par- 
allelized for three-dimensional applications. We imple- 
mented the transport of UV plane waves along the x-axis 
of a cartesian box (or along the z-axis of a cylindrical co- 
ordinate grid) in addition to the original radial transfer 
of photons from a point source centered in a spherical 
polar coordinate mesh. The parallel reactive flow, addi- 
tions to the radiative transfer, and improvements to our 
subcycling scheme are detailed below. 

2.1. Parallelization Scheme 

We refer the reader to WN06 for a description of our 
sequential multistep radiation algorithm. ZEUS-MP de- 
composes the simulation domain into subunits called 
tiles, with each tile assigne d to a single process or for the 
duration of a calculation ijHaves et al. I I2OO60 . Global 
solution concurrency is maintained by message-passing 
between processors through calls to the Message Passing 
Interface (MPI) library. In general, the computational 
domain can be partitioned along all three axes but in 
photoionization calculations we limit the decomposition 
to two axes only. In cartesian coordinates we subdivide 
the volume along the y and z directions only; in cylin- 
drical coordinates the decomposition is along the r and 
phi axes, and in spherical grids we tile only along theta 
and phi. Hence, a tile in a cartesian grid is a rectangu- 
lar box whose length in the x direction spans the entire 
grid. In cylindrical coordinates, tiles are wedges in r and 
phi lengthwise in z, and in spherical coordinates they 
are slices in theta and phi radiating outward from the 
coordinate center to the outer radial boundary. 

Domain decomposition is implemented along only two 
coordinates in ZEUS-MP to enforce load balancing of the 
radiative transfer and to avoid the need to communicate 
photon fluxes between adjacent tiles. The key to effective 
load balancing is requiring that the ionization front be 
present in all the tiles at once. Subdivision along all three 
axes for a single point source or plane wave would result 
in the few tiles hosting the front shouldering a dispropor- 
tionate fraction of the chemistry subcycling, leaving the 
other tiles relatively idle. The tradeoff in this scheme is 
the limit to the parallelization that can be applied to a 
given problem. 

The reaction network is evolved plane by plane within 
a tile. Updates to the reaction network can be made 
one plane of zones at a time because the chemistry in 
any given zone over Atchem does not depend on its near- 
est neighbors. We advance the network one plane at a 
time so the reaction coefficients can be stored in two- 
dimensional arrays. Up to thirty rate coefficients are 
present in the reaction network and as many as 21 heat- 
ing and cooling coefRcients comprise the isochoric update 



to the gas energy. Any network update order requiring 
the coefficient arrays to be three-dimensional would in- 
flate the memory required for a calculation by a factor 
of ten beyond what is required just for hydrodynamics. 

We prevent the solution for the entire grid from being 
evolved further than the shortest evolution time in any 
of its tiles by performing an MPI alLreduce operation to 
extract the global minimum of Atheat- The smaller of the 
grid minima of Atheat and the Courant time is adopted 
for the hydrodynamical update of gas energies, densities, 
and velocities. To further preserve solution coherence, 
mpi_barrier calls are stationed at the end of the reaction 
network module to ensure subcycling is complete in all 
the tiles before any of them advance to the next operator- 
split step in the hydrodynamics. This procedure is nec- 
essary because of the implementation of asynchronous 
message passing in ZEUS-MP, in which code execution 
in some tiles can proceed while communication between 
other tiles is still in progress. 

2.2. Radiative Transfer: Cartesian and Cylindrical 
Coordinates 

We compute the attenuation of the ionizing flux along 
rays parallel to the grid lines. In cartesian coordinates, 
the rays are parallel to the x-axis. In cylindrical coor- 
dinates, the rays are parallel to the z-axis. Generalizing 
the radiative rate coefficients k^ad evaluated in WN06 to 
planar radiation flows is straightforward. Recalling that 
the number of photons removed from a zone per second 
i'n-ioniz) is related to the radiative rate coefficient (krad) 
by eq. 18 of WN06 

(1) 
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and that photon conservation demands that 
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where F and A are the flux (in erg/sec/cm^) and area 
over the given face of the zone, we have that 
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where hph is the emission rate of all photons at a given 
energy. Ax is the zone length along the x-axis, and x 
is the mean free path of the photons in the neutral gas. 
This derivation holds with axial plane wave transport in 
cylindrical geometries because the metric terms in the 
area factors drop out. We include an option for atten- 
uating the plane wave flux by r~^ to approximate the 
geometric dilution of a point source. 

We tested the radiative transfer by computing the posi- 
tion of a one-dimensional planar ionization front in a uni- 
form static medium of hydrogen, with n^f — 1000 cm~'^ 
and T = 72 K. The grid has 1000 zones in the x direction 
and a length of 1.0 pc. The photon flux along the x-axis 
incident to the left-hand face is Jq = 4.0 x 10^^ cm~^ 
s~^, and we applied a constant case B recombination co- 
efficient aB = 2.015 X lO^^'^ cm"^ s^^. By equating the 
photon flux at the front to the incoming flux of neutral 
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Fig. 1. — Approach of the static planar ionization front to xgtr, 
numerical vs analytical. 



particles we can compute the velocity and position of the 
front in the usual manner: 
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where n^j is the total number of ionizing photons reach- 
ing the front per unit time and A is the area of the zone 
face. The position of the front easily follows from simple 
integration: 
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We show the evolution of the front in the static medium 
in Fig[l] for the parameters of this test, y^str — 0.644 pc. 
The code agrees with the analytical prediction to within 
5% at all times. A parallelized version of this test was 
performed in three dimensions to verify that the advance 
of the front in each tile was identical. 

2.3. Advanced Subcycling 

As detailed in WN06, in one-dimensional problems the 
primordial gas reaction network in ZEUS-MP was cycled 
over consecutive chemical timesteps 
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until a hydro timestep was covered 
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This scheme proved to be robust and versatile, accommo- 
dating regions experiencing either photoionization or re- 
combination. However, timescale analysis revealed that 
cells close to the source without any initial electron frac- 
tion required as many as 70 - 80 hydrodynamical up- 
dates before coming to ionization equilibrium upon ar- 
rival of the front. Network subcycling between updates 
also tended to be extreme in such zones, with several 



thousand iterations required to cross At^ydro at times, 
far more than are required for accuracy. The problem is 
ameliorated at larger radii because heating times increase 
further away from the source. Also, advection transports 
small electron fractions into zones slightly ahead of the 
front, preventing initially minute Atchem- 

Practical Tfront transport in three dimensions required 
two improvements to our integration scheme. When 
the front encounters overdensities, advection of electrons 
into those zones is halted, resulting in paralyzingly short 
chemical timesteps until the front exits the perturbation. 
We eliminated this problem by modifying Atchem- 
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The correction term ensures that At^ 



is not too small 



when photoionizations commence in the zone and has lit- 
tle effect as the zone approaches equilibrium and Atchem 
increases. The algorithm now advances ionization fronts 
through density jumps without difficulty. 

We also experimented with more agressive subcycling 
by requiring that a zone come to ionization equilibrium 
in no more than 7-8 hydrodynamical updates by altering 

At fly dro- 



At 



hydro 
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^ht/cl 



(9) 



We obtain speedups of up to 15 in transporting R-type 
fronts with this time step prescription, but there is little 
improvement with D-type fronts because in primordial 
environments Atcour is usually the smaller of the two 
timescales in eq. |9l Two difficulties arose in this method. 
First, the longer time step can cause overshoot in the 
cooling as the zone reaches full ionization, occasionally 
removing all the gas energy in a zone in a single update. 
We resolved this issue by reverting to the original time 
step if igas < in a zone, which removed the anomalous 
cooling with little sacrifice in performance. 

Unfortunately, the coarser time steps can lead to inac- 
curacies in the reaction network, causing the propagation 
of R-type fronts to diverge from theory by as much as 10 
- 20%. In many practical applications this is not a se- 
rious error, so the performance gain may outweigh the 
loss in accuracy. While we retain both subcycle options 
in the code, we applied only the improved Atchem to 
these studies to properly capture peturbations in R-type 
fronts that later grow into instabilities after becoming 
D-type. 

3. D-TYPE INSTABILITIES IN TWO DIMENSIONS 

We now examine the growth of unstable modes in a 
D-type ionization front that is expanding in a radially 
symmetric gas profile randomly seeded with minute fluc- 
tuations in density on a polar coordinate grid. We con- 
sider the development of instabilities in a uniform core 
enclosed by an r~^ envelope 



Uc if r < Tc 

ndr/rc)^'^ if r > Tc 



assuming that the shocked gas is radiatively cooled by 
metal emission lines. This problem was first studied by 
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iGarcia-Segura fc Francol (|1996f ) (hereafter GSF) and we 
compare their results below. 

We augmented the isochoric gas energy update in 
ZEUS- MP (eg. 9 of WN06) with cooling curves 
from Dalgarn o fc McCravl ()1972f ) in which any elemen- 
tal abundances can be specified. The Dalgarno & Mc- 
Cray curves include both electron and neutral hydrogen 
collisional excitation of H, He, and neutral and singly- 
ionized C, O, N, Fe, Si, and S. The energy equation also 
has cooling due to electron collisional ionization of hydro- 
gen, H and He recombinational cooling, bremsstrahlung 
cooling, and Compton cooling by free electrons off the 
cosmic microwave background (for high-redshift applica- 
tions). These latter two processes play a negligible role 
in our simulations. In general, the Dalgarno fc McCray 
cooling rates in a given cell can be constructed with any 
electron or neutral H number density derived from the 
reaction network. For consistency with GSF, we evaluate 
them assuming a fixed electron fraction of 0.01 and solar 
metallicity in the calculations in this paper, regardless of 
the true electron abundances on the grid. 

This treatment approximates radiative cooling in 
shocked neutral gas well but underestimates metal line 
cooling, H-a cooling, and He-a cooling in ionized gas 
because the electron fraction is much higher than 0.01 
there. At postfront gas temperatures above ~ 10^ K, 
H-a and He-a cooling rates per atom can exceed those 
of the metals by four orders of magnitude. Given that 
neutral H and metal fractions in the ionized gas at solar 
metallicities are roughly equal (~ 10~^), H-a and He-a 
cooling dominates the total Dalgarno fc McCray curve 
inside the H II region. However, in the pure H simula- 
tions of this study we incur overestimates of at most 10% 
in temperature in the ionized gas because it is H recom- 
binational cooling that really controls postfront temper- 
atures, as we demonstrate below. Cutoff temperatures 
were also incorporated in the runs, defined to be the gas 
temperature below which cooling is not applied to the gas 
energy density. For the two runs discussed in this sec- 
tion we adopted cutoff and background of temperatures 
of 100 K. 

A point source of ionizing UV flux was centered on 
the polar coordinate grid with 200 radial zones and 180 
zones in the theta direction. The inner and outer ra- 
dial boundaries were 0.01 pc and 2.0 pc with reflecting 
and outflow boundary conditions, respectively. The in- 
ner and outer boundaries in the theta direction were 7r/4 
and 37r/4, respectively, with reflecting boundary condi- 
tions. We adopted this angle range to avoid the minute 
Courant times associated with zones near the coordinate 
poles due to their small angular dimensions. At the res- 
olution of this simulation we obtain minimum Courant 
times ten times greater over this range than at the poles. 
Our choice of reflecting boundary conditions could ex- 
clude the growth of long- waveleng t h mod es, but the lin- 
ear stability analysis of iGiulianil (|1979[ ) predicts that 
short-wavelength modes will initially dominate and then 
later consolidate into larger structures. 

The density profile in both runs was hydrogen only, 
with a flat central core of 1.0 x 10'' cm~^ out to r = 0.2 
pc followed by an r~^ dropoff. At the outer boundary the 
density falls to ~ 100 cm~'^. The initial gas temperature 
was set to 100 K. We adopted a photon emission rate riph 



of 1.0 X 10'*^ s^' with energy 15.2 eV, with an energy per 
ionization er of 1.6 eV in order to maintain the ionized 
gas temperature at 1.0 x 10'* K. Theory predicts that if 
the initial Stromgren radius falls in the density gradient 
the I-front will re vert from D-type to R-type and over- 
run the envelope (jFranco et al. lfT990t ). For the central 
density and emission rate of these runs, the Stromgren 
radius 

/ o • \ 1/3 

is 0.065 pc, entirely within the central core. The transi- 
tion to R-type is thus delayed. 

In the first run we do not perform radiative transfer 
but instead repeat the simplified treatment applied in 
the GSF models. Every hydrodynamical time step the 
equilibrium position of the I-front is computed along a 
given radial ray by equating the sum of the case B re- 
combinations along the ray to the emission rate of the 
central source: 

/ A-Kr^jieripaBdr = fiph (11) 
Jo 

The gas temperature is then set to 1.0 x 10* K along 
the line of sight out to r/. After repeating this proce- 
dure for every value of theta the model is evolved for a 
Courant time. The position of the front thus changes 
each time step as gas is rearranged on the grid by the 
high-temperature fiow. We refer to this approach in the 
following as "ionization equilibrium" . In practice, we 
found it necessary to limit the cooling in a shocked zone 
to be no greater than 90% of the energy remaining in the 
zone above the temperature cutoff of the cooling curve. 
Otherwise, the rather coarse hydrodynamical timesteps 
can allow explicit updates to the energy equation to com- 
pletely remove all energy within the cell. No R-type 
fronts ever arise in this approximation since the fronts 
by construction are always assumed to have reached the 
Stromgren at any given time step. 

In the second run we transport the front across the 
grid using radiative transfer scheme described in Sec. 
2.2, coupled to self-consistent photoionization and heat- 
ing as_described_Jn_Sec^_^2^__(see also WN06). In 
the iGarcia-Segura fc Francol (|1996l ) simulations, density 
fluctuations with randomly-distributed amplitudes of at 
most 1% were imposed on very zone to seed the forma- 
tion of shock instabilities. We also did this, holding the 
gas energy density constant to prevent pressure fluctua- 
tions from dispersing the perturbations. However, in our 
model we only imposed these variations beyond radii of 
0.125 pc to prevent the onset of shadowing instabilities 
in the early R-type front before dynamical instabilities 
in the D-type front can form. We compare the evolution 
of both models in Fig[2l 

3.1. Ionization Equilibrium vs Radiation 
Hydrodynamics 

Both models confirm the prediction of IGiulianil (119791 ) 

that the modes which grow fastest are those with the 
largest wavelengths. Unstable modes begin to disrupt 
the shock at 9.5 kyr in the second run, with 13 spikes vis- 
ible at 28.5 kyr (panel (b)) but only 5 remaining at 148.4 
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Fig. 2. — Dynamical instability in a D-type I-front propagating out of an density gradient with metal line cooling. Panels (a), (c), 
and (e): number densities at 28.5 kyr, 76.1 kyr, and 148.4 kyr for a front evolved assuming simple ionization equilibrium. Panels (b), (d), 
and (f): density evolution of an I-front transported across the grid with full radiative transfer at 28.5 kyr, 76.1 kyr, and 148.4 kyr. 



kyr (panel(f)). Vorticity in the flow is evident where the 
bases of the spikes join the flared regions. The instabih- 
ties remain D-type: they are always bounded by a dense 
shocked gas layer and their ionized interiors have time 
to laterally expand as they elongate. The perturbation 
amplitudes never saturate in either model but continue 
to grow well into the nonlinear phase, from 0.1 pc at 28.5 
kyr to 0.5 - 1.0 pc at 148.4 kyr in the second run. Direct 
quantification of their growth rates is problematic be- 
cause individual structures coalesce over time, but their 
rapid elongation at later times is due to their descent 
down the density gradient. 



There are noticeable differences in structure of the 
fronts at early and intermediate times. While the ionized 
fingers are numerous and narrow in both simulations at 
early times (panels (a) and (b) of Fig [21), their morpholo- 
gies diverge after a few tens of kyr. Ionized spikes push 
past each other at random times through the shock in the 
second run, flaring sideways earlier than in the first run 
because they are not as collimated by their neighbors. 
In contrast, the spikes advance nearly in step with one 
another in the first model, remaining cylindrical out to 
much greater distances because of the parallel advance of 
their neighbors. There are also less of them: 7 spikes at 
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Fig. 3. — One-dimensional I-front positions in the ionization 
equilibrium approximation and with radiative transfer 



0.3 pc in comparison to 13 in the second simulation. The 
differences between the two models are most pronounced 
at early and intermediate times, becoming less distinct 
as they enter highly nonlinear regimes. 

The departures are chiefly due to how the two numeri- 
cal schemes capture the interaction of the front with the 
shocked neutral shell. The topology of the instability 
is sensitive to how the inner shell is photoevaporated. 
Ionization equilibrium abruptly assigns temperatures of 
1 X 10^ K to the fluid elements on the inner surface 
that are capable of being ionized. Had they instead been 
incrementally heated and ablated, as in the radiation hy- 
drodynamical scheme, their subsequent expansion would 
have created a larger number of cracks in the shell dis- 
tributed over time, yielding the larger number of ionized 
fingers observed in the second run. Although costlier, 
full radiation transport is necessary to model the actual 
deformation of the shocked shell. Initial disparities in the 
deformation of the shell by the two algorithms lead to the 
much greater differences in morphology at intermediate 
times evident in panels (c) and (d) of Fig [21 

In one-dimensional calculations the two algorithms 
predict I-front positions to within a few percent of each 
other at early times but diverge by more than 10% at 
later times, as shown in Fig [31 This agreement holds only 
so long as the front is D-type. Ionization equilibrium ap- 
proximations would not fare as well if the front became 
D-type in the r~^ gradient or in other stratified media. 
It would fail altogether if the front executes transitions 
from D-type to R-type or vice versa, as is common in 
cosmological H II regions. We have found that updates 
to velocities should be performed on photoheating time 
scales rather than the Courant time to capture the true 
acceleration of fluid elements by the front in power-law 
density gradients (WN06). 

4. THREE-DIMENSIONAL D-TYPE INSTABILITIES 



We extended the radiative transfer calculation in sec- 
tion 3 to three dimensions using a spherical polar coor- 
dinate grid with 200 radial zones and 180 zones in the 
theta and phi directions. The inner and outer radial 
boundaries were again 0.01 pc and 2.0 pc with reflecting 
and outflow boundary conditions, respectively. The an- 
gle boundaries were again chosen to be 7r/4 and 37r/4 in 
phi and theta with reflecting boundary conditions. We 
apply the same randomly perturbed density profile as 
before to the grid, with the energy per ionization ep set 
to 1.6 eV to ensure a postfront gas temperature of 1.0 x 
10* K. 

Our first run (S31) had cutoff and background tem- 
peratures of 1000 K and 10 K, respectively, while our 
second run (S22) had cutoff and background tempera- 
tures of 100 K. Metal lines cannot fully cool shocked gas 
in the first case so the shell is thinner than in an adia- 
batic flow but not as dense (or cold) as when the metals 
are free to cool the gas down to t h e bac kground tem- 
perature. iGarcia-Segura fc Francol {1996) nevertheless 
found instabilities to readily form in these circumstances 
in their two-dimensional simulations. The lower cutoff in 
the second run is intended to produce a relatively cold, 
dense shell that will break up more violently. In both 
of these models fiph was again chosen to be 1.0 x 10'*^ 
s~^ so that the Stromgren sphere would reside within the 
core. 

4.1. S31 Instability 

In Fig ([4]) we show (f> — Tr/2 slices of the evolution of 
the ionization front in the S31 run at 29.5 kyr, 70 kyr, 
and 100 kyr (the total time for which the front was fol- 
lowed in the 1996 experiments). In panels (a) and (b) 
the perturbations are seen to be dominated by short- 
wavelength modes at early times that later merge into 
the longer- wavelength modes in panels (c) and (d), in 
accord with linear stability analysis. The initial wave- 
lengths are of the order of the shell thickness. Unable 
to completely cool, the shocked shell becomes thicker as 
neutral gas accumulates onto it, acquiring turbulent mo- 
tions as it evolves. The instabilities persist to 200 kyr in 
these three-dimensional calculations as seen in panels (e) 
and (f) but saturate without the prominent growth ex- 
hibited in the GSF two-dimensional models. The turbu- 
lence in the shell appears to prevent the unstable modes 
from breaking out. The maximum density attained by 
the clumps in the course of the simulation is 4.99 x 10^ 
cm~"^; pronounced clumping is visible at 177 kyr in panel 

(f). 

4.2. S22 Instability 

The ionization front breaks through the shock and es- 
capes much further down the stratified envelope in the 
S22 run, as dramatically illustrated in Fig We show 
(j) ~ tt/2 slices of the front for this model at 20.8 kyr, 67.9 
kyr, and 124.5 kyr. Radiation fiares outward through 
the cracks in this more violently fragmenting, colder, 
denser shell, and the dominant growth modes over time 
are again those with the greatest wavelengths. This mo- 
tion is visible in the large dark neutral patches cutting 
across some of the ionized flares in temperature panel (e) 
of Fig ([5]) . Comparison of panels (c) and (e) demonstrate 
the transience of these features. Newly ionized gas in the 



Fig. 4. — Dynamical instability in a D-type front propagating outward along an r ^ density gradient with metal line cooling, run S31. 
Panels (a), (c), and (e): temperature evolution at 20.8 kyr, 96.2 kyr, and 177.4 kyr. Panels (b), (d), and (f): density evolution at 20.8 kyr, 
96.2 kyr, and 177.4. 



outer regions of the flares expand laterally as the insta- 
bilities propagate, creating dark, rarefied zones in the 
corresponding density panels. The unstable modes ad- 
vance as D-type, given their time scales of expansion as 
well as because they remain bounded by shocked gas that 
is visible as a thin, lighter hued density layer in the outer 
edges of the flares. The maximum density of the gas in 
this run is 2.08 x 10^ cm~'^, greater than in the previous 
run due to the higher degree of postshock cooling. The 
thin radial segments of cool zones visible in temperature 
panels (c) and (e) are cold, dense shell fragments of the 
shocked shell not yet ionized by the front. These frag- 



ments can be seen to originate in the cold shell because 
they extend from the base of the instabilities and radi- 
ate outward. They are overtaken on either side by hot 
ionized outflows, with the neutral parcels of gas closest 
to the central source eventually becoming ionized them- 
selves. Their large densities stand out in relief against the 
surrounding dimmer regions in panels (d) and (f). The 
cooling of the fragments down to the cutoff temperature 
is due their high densities rather than some instability in 
the explicit updates to the energy equation. Our integra- 
tion scheme prevents any energy loss greater than 10% 
over the time step on which hydrodynamic updates are 



Fig. 5. — D-type ionization front instability in an r ^ density gradient with metal line cooling, run S22. Panels (a), (c), and (e): 
temperature evolution at 20.8 kyr, 67.9 kyr, and 124.5 kyr. Panels (b), (d), and (f): density evolution at 20.8 kyr, 67.9 kyr, and 124.5 kyr. 



performed. Large clumps of gas can be seen to persist 
to late times as they are blown outward by the front in 
panel (f). 

There are clear qualitative differences between the un- 
stable flow in the two-dimensional radiation hydrody- 
namical model in the last section and this run that can 
only be attributed to their dimensionality, given that 
both otherwise employ the same physics. First, trans- 
verse flow of gas into the plane of the image from above 
and below is evident in some of the dark patches in the 
outer regions of the ionized fingers in panels (c) and (e) . 
These are not cold dense shell fragments (as the density 



panels attest) but instead are neutral gas crowded into 
the plane of the image by adjacent instabilities. Sec- 
ond, the shapes of the fingers differ in the two runs: the 
structures in the three-dimensional model are both more 
narrow and plentiful than in the two-dimensional run. 
Less flaring occurs in the three-dimensional fingers be- 
cause they are partially coUimated by neighboring flows 
that are not possible in the two-dimensional models. The 
geometry of cracks in the shell at the base of the insta- 
bilities may also be altered by the additional degree of 
freedom in the flow. 
The S31 and S22 data confirm that Lfront breakout 
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Fig. 6. — Hydrodynamical profiles for ionization fronts in a spherically-symmetric r ^ density stratification. Solid lines are I-front profiles 
with gas cooling truncated at 1000 K while dotted lines denote profiles in which gas cooling is halted below 100 K. Left panel: Rg < Vc- 
1000 K profiles are at 0.48 Myr, 0.82 Myr, and 1.37 Myr while the 100 K profiles are at 0.7 Myr, 1.01 Myr, and 1.26 Myr. Right panel: 
Rs > rc. 1000 K profiles are at 36 kyr, 40 kyr, and 44 Kyr while the 100 K profiles are at 32 kyr, 38 kyr, and 42 kyr. 



from the envelope is determined by both the degree of 
fragmentation of the sheU (which is governed by the ef- 
ficiency of coohng in the shocked gas) and the pressure 
ratio of the ionized and background gas. At early stages 
of the instability there are clear differences in the initial 
scales of the perturbations that we suspect, but have not 
confirmed, are due the thickness of the shocked gas shell 
at the time of initial fragmentation. The shell thicknesses 
vary from five to ten zones at early times, within the nu- 
merical resolution of our mesh. The initial scales of these 
deformities in part differentiate the later topologies of the 
S31 and S22 instabilities. 

We further note that the absence of runaway instabil- 
ities in our S31 model may in part be due to R5 being 
less than ic- In one dimension, somewhat surprisingly 
we find that the density envelope actually confines the 
front behind the shock at distances far beyond the prob- 
lem boundaries in our calculations, in contrast to the 
relatively quick breakout that would occur had the ini- 
tial Stromgren radius fallen within the gradient. It is 
therefore likely that the perturbed density profile would 
prevent further growth of the instabilities until they ap- 
proach the breakout radii found in our one-dimensional 
calculations discussed below. Accurate models of I-front 
transitions are therefore key to capturing the correct 
growth of instabilities and demand full radiation hydro- 
dynamical treatment of the ionized flows. 

If the Stromgren radius lies beyond the flat central core 
in the one-dimensional density profile of eq. [3] then its 
value is given by 



Rlu — Rs 
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(12) 



where Rs is the Stromgren radius for the front in a uni- 
form medium. If > r^, the I-front executes the classical 
transition to D-type but the n soon reverts t o R-ty pe as 
it descends the gradient. The [Franco et aC] (|1990l ) anal- 



yses do not address the evolution of the front when it is 
initially confined to the core (R5 < rc), as in these mod- 
els. To determine the outcome of the front in this regime 
we performed two sets of one-dimensional calculations 
in the density profile used earlier but without perturba- 
tions. We extended the outer boundary of the mesh to 
20 pc with 2000 radial zones to preserve the spatial res- 
olution. One set had a central photon rate fiph = 1.0 
X 10** s^^ as before while hph — 5.0 x 10*^ in the 
other set. This larger rate sets R^ — 0.26 pc, outside the 
central plateau. One calculation in each set had a cut- 
off temperature of 1000 K while the other had a cutoff 
of 100 K; the background temperature in each was 100 
K. The evolution of the front was followed for 1.3 Myr 
in the first set and for 45 kyr in the second set (much 
less time is required when there is rapid escape down the 
gradient). Flow profiles for all four calculations appear 
in Fig [6l 

Even out to 20 pc, ten times the distance to which our 
earlier models were followed, the density and ionization 
fraction profiles in the left panel of Fig [6] reveal that the 
I-front remains trapped by the shock at the lower pho- 
ton rates. The 100 K shell is only half as thick as the 
1000 K shell but has ten times its density. The higher 
recombination rate in the cooler shell offsets its smaller 
width to confine the front. Note that radiative cooling 
in both cases lowers the shell to the cutoff temperature. 
When we extend the grid boundary to 200 pc at the 
same spatial resolution (thereby ensuring that the front 
remains on the grid fo r the 10 Myr lif etime of a star with 
these emission rates (|Schaereij|2002[ )). we find that the 
front finally breaks through the 1000 K shell at 90 pc 
and through the 100 K shell at 64 pc. In these calcula- 
tions the shock gains strength as it descends the density 
gradient. Line cooling still collapses gas into a relatively 
thin layer at the base of the shell but cannot cool all the 
material that accumulates as the shock accelerates so the 
shell heats and thickens. When the shell's outer layers 
reach 10,000 K, coUisional ionizations assist the I-front 
to break through the shock. 
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Fig. 7. — Run S31 with R^^ > Vc- Panels (a), (c), and (e): temperature evolution at 9.4, 13.2, and 15.1 kyr. Panels (b), (d), and (f): 
density evolution at 9.4, 13.2, and 15.1 kyr. 



On the other hand, the right panel demonstrates that if 
Rs = 0.26 pc > Tc, the f ront breaks free of th e shock just 
past Rg as predicted by lFranco et al. I ()199Clf ). Reverting 
to R-type, the front quickly exits the spherical cloud, 
leaving behind a steepening ionized core shock at the 
site of detachment as seen in the density panel. Little 
gas dynamical evolution occurs over the time in which 
the front crosses the grid. The higher density of the 
cooler shell slightly delays I-front breakout (by ~ 2 - 
3 kyr). This implied that instability growth would be 
particularly explosive in cases where Rs lies beyond the 
core radius. 



4.3. S31 Breakout Instability 

To test this conjecture we devised a new S31 run with 
nph = 5.0 X 10"*^ s^^. We imposed random variations on 
the density as before but beyond radii of 0.29 pc in order 
to allow the I-front to become fully D-type at R^^ — 0.26 
pc before encountering the perturbations. Temperature 
and density images for (j) = tt/2 in this calculation are 
shown in Fig ([7]) for 9.4, 13.2, and 15.1 kyr. Radiation 
breakout through cracks in the shock rapidly ensues and 
the front erupts through the envelope in spikes which 
are R-type, as evidenced by the short escape times and 
the lack of hydrodynamic response in the ionized fingers. 
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The slightly jagged edges in some of the spikes are arti- 
facts of the spherical polar coordinate mapping to carte- 
sian coordinates. Fragmentation of the shell is visible at 
early times as bright dots in the shell in the density pan- 
els but it soon disappears as the clumps are photoionized 
and then dispersed on time scales setby the sound speed 
in the postfront gas. The maximum density achieved in 
this simulation was 5.47 x 10'' cm~'^. The spikes remain 
fairly narrow because they elongate rapidly in compari- 
son to hydrodynamical time scales, but they do begin to 
drive weak transverse shocks at their base as the ionized 
gas begins to laterally expand. These unstable modes do 
not fall in any of the categories discussed so far because 
they evolve as fast R-type structures. 

5. R-TYPE SHADOW INSTABILITIES 

In this test, a plane-parallel R-type ionization front 
enters the yz-face of a rectangular box with a uniform 
gas density except for a spherical bump slightly offset 
from the face of entry. The density in the bump varies 
linearly in radius from the ambient value at its center to 
either 50% below or above this value at its surface. Our 
simulation volume has a length of 1.0 pc along the x- 
axis and 0.25 pc along the y and z-axes, with 1000, 250, 
and 250 zones in the x, y, and z directions respectively. 
The gas number density in the box is 1000 cm~'^ and is 
hydrogen only. The temperature of the gas was set to 
72 K, yielding a sound speed Cg of 1 km/s. The fixed 
energy er per photoionization was 0.8 eV to yield initial 
postfront temperatures of 10^. The incident photon flux 
along the x-axis at the entry face was 3.0 x 10^^ cm^'^ s^^ 
to approximate that of an O star at its Stromgren radius. 
The radius of the perturbation was 0.0125 pc and was 
positioned 0.05 pc along the x-axis and centered in the 
yz-plane. Reflecting and outflow boundary conditions 
were applied at the pc and 1.0 pc faces, respectively, 
with reflecting boundaries along the other four faces. 

Our physical parameters for this test are similar to 
those of the two-dimensional simulations of IWilliami 
(|1999f ). except for cooling. In those models near- 
isothermality was enforced by setting 7 = 1.1 and the 
internal energy of the ionized gas to be constant to guar- 
antee Cs ~ 10 km/sec. In contrast, we applied the same 
cooling algorithm described earlier with a cutoff temper- 
ature of 1000 K. The front reaches a Stromgren distance 
of 0.32 pc at approximately 1000 yr: 



where Fq is the photon flux (cm^^ s^^), n is the number 
density of neutral hydrogen and is the case B re- 
combination coefflcient. The simulation volume was par- 
titioned into 25 uniform tiles, with five divisions along 
both the y-axis and z-axis and each tile spanning the 
length of the volume in the x direction. 

We show xy slices at z = of the temperature and den- 
sity evolution of the shadowing instability in Figs [8] and 
[9] due to an underdense perturbation at 1.9 kyr, 4.9 kyr, 
and 12.6 kyr. At 1.9 kyr the front has transitioned to 
D-type and the original corrugation in the R-type front 
is well into the nonlinear growth phase. The expand- 
ing photoevaporated perturbation is faintly visible on the 
left of the figure. The twin peaks in the instability are 
xy slices through a ringed jet that is due to the radial 



density profile of the spherical bump: the R-type front 
preferentially advances along the lines of sight parallel to 
the X-axis that cut the underdense regions close to the 
surface of the bump. Along the line of sight piercing the 
bump through its center, the front advances at nearly 
the same rate as in the unperturbed medium because 
the densities along this path through the bump are close 
to those in the surrounding gas. 

At later times cooling in the shell fractures the insta- 
bility into multiple jets with highly nonlinear growth and 
morphologies reminiscent of the density structures visi- 
ble in the two-dimensional numerical experiments. No 
temperature images of the 1999 work are available so it 
is not possible to directly compare the distribution of 
ionized and shocked gas within the jets or verify the ex- 
istence of the ionized extensions into the shock that are 
present in our models. The unperturbed regions of the 
front snowplow gas in to a somewhat broader shell than 
in the iWilliami ()1999D models because the gas can only 
cool to 1000 K. The length of the jet in our model (0.14 
pc) is somewhat shorter than in the Williams runs due 
in part to the lower ionizing flux. However, intrinsic dif- 
ferences between two and three-dimensional flows may 
also come into play, as the additional degree of freedom 
opens transverse flow channels that may quench the in- 
stability. Indeed, some vorticity is evident in the flow at 
later times, a possible signature that instability growth 
is coupling to turbulent motion in the gas and being re- 
duced. The densities in the jets range from 19.0 cm~'^ 
to 3.72 X lO'' cm~^, a factor of 4 or 5 lower than in 
the 1999 experiments. The greater maxima in the earlier 
work in part arise from the higher density jumps permit- 
ted across the isothermal shocks in those models as well 
as from the lower temperatures to which the postshock 
gas can cool. As the central jet advances it widens as 
ionized gas in its interior expands perpendiculary to the 
flow. 

Two features also worth mention are the density knots 
in the shock located above and below the midplane of the 
jets at their base at 1.9 kyr. As the front advances, radi- 
ation escapes past the knots and forms jets that migrate 
across the face of the shock as they elongate. Shortly 
thereafter, smaller disturbances in the shock materialize 
further from the base of the central jet, possibly triggered 
by acoustic waves along the shock or perhaps merely nu- 
merical artifacts. These density fluctuations quickly cool 
and create cracks in the shock through which ionizing 
photons penetrate, forming smaller subsidiary jets that 
later merge with the widening central jet in complicated 
hydrodynamical flows. The original knots are present in 
the two-dimensional runs and are attributed by Williams 
to a type of zero- wavelength odd-even numerical instabil- 
ity (zero-wavelength in the sense that it would be mani- 
fest at all resolutions). The phenomena apparently arise 
because the I-front itself is not resolved by the grid. Just 
as postshock oscillations can come to dominate shocked 
flow not broadened by a few zones with artificial viscos- 
ity in lower-order hydrodynamical schemes, under certain 
circumstances sharp I-fronts can preferentially advance 
in an oscillatory fashion across adjacent zones. Williams 
removed this effect in his simulations by setting the pho- 
toionization cross section to be a function of zone width 
(in direct analogy to artificial viscosity), hence broad- 
ening the front. Practical applications employing mul- 
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Fig. 8. — Temperature evolution of an I-front shadow instability due to an underdense spherical underdensity; (a) 1.9 kyr, (b) 4.9 kyr, 
and (c) 12.6 kyr. 



tifrequency radiative transfer (for Population III stellar 
and miniquasar hard spectra) may obviate this feature 
because the spread in photon mean free paths would nat- 
urally widen I-fronts in many environments without re- 
course to artificial means. 

6. CONCLUSION 

For the past ten years, the general paradigm of ioniza- 
tion front instabilities in the ISM has been one of sword- 
like structures with unrestrained amplitude growth in a 
wide variety of environments. Our new models temper 
this view with saturated modes that likely drive turbu- 
lent flows in some instances. The morphology of the S31 
instability in particular suggests that some modes may 
efficiently drive turbulence in the shocked shell, with pos- 
sible implications for time scales of subsequent star for- 



mation in the region. The numerous short-wavelength 
ionized spikes in the shock at early times later gather 
into fewer long- wavelength structures, so energy is trans- 
ported at first from smaller scales to larger ones by the 
flow. Later, these large scale flows crumple into ever 
smaller scales as the saturated instabilities devolve into 
shear motions and develop the vorticity evident in the 
Kelvin-Helmholtz features appearing in panel (f) of Fig 
[H We expect this process to be more prominent in actual 
molecular clouds because the r~^ proflle of a core trun- 
cates in relatively high intra-cloud densities that would 
further flatten the instabilities. While our current mesh 
can resolve the shell breakup and ionized fingers well, 
these features would have to be evolved for longer times 
on finer grids able to capture turbulent cascades to deter- 
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Fig. 9. — Density evolution of an I-front shadow instability due to the spherical underdensity: (a) 1.9 kyr, (b) 4.9 kyr, and (c) 12.6 kyr. 



mine the extent to which turbulence arises in the cloud. 
Nevertheless, the possibility that ionization front insta- 
bilities together with winds from massive stars and su- 
pernova explosions may provide an important source of 
turbulent support in star forming clouds in the galaxy 
today is an intriguing one. 

We find that the dimensionality of our models distin- 
guishes them from previous work more than their im- 
proved physics. The structures emerging in our simu- 
lations bear resemblance to many objects visible in the 
ISM today. The globules and pillars in the S31 run are 
reminiscent of those seen in the Eagle nebula while the 
long wispy extensions in the S22 run are comparable to 
those in many planetary nebula such as NGC 6751. The 
clumping that is prevalent in the galaxy is also evident 
in our models. In reality, ionization front instabilities 



compete with wind-blown structures and blast waves to 
produce the panoply of features in the local interstel- 
lar medium. How all three act in concert to create the 
circumstellar environments of massive stars and gamma- 
ray bursts is an important question to be addressed in 
future work. Likewise, how H II region instabilities inter- 
act with the magnetic fields threading molecular clouds 
remains poorly understood. We expect that unstable 
modes will not be seriously suppresed on parsec scales 
because of the relative magnitudes of ionized gas and 
magnetic field pressures, but their evolution could be al- 
tered on larger scales as they approach pressure equilib- 
rium with their surroundings. 

A question of relevance to cosmological reionization is 
whether I-front instabilities arise in the primordial clouds 
of primeval massive stars and protogalaxies. In pristine 
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gas shocks must resort to atomic hydrogen line cooling 
and Compton cooling at high redshifts, which are far 
less efficient processes than metal line emission. If the 
shocked flows are non-radiating, no Vishniac instabili- 
ties can develop or be amplified by radiation so the fronts 
would be stable. On the other hand, 1-front shocks ac- 
celerate if they emerge from cosmological minihalos with 
density gradients steeper than r~^. Their exit may thus 
be subject to Rayleigh- Taylor instabilities that ionizing 
radiation could blossom into much larger structures. We 
further note that hard Population III or miniquasar UV 
and x-ray spectra can significantly broaden comological 
I- fronts, and it is unclear what impact the layer of partial 
ionization will have on instability formation. These ques- 
tions together with how 1-fronts sculpt gas clouds in the 
early intergalactic medium will be the focus of studies to 
come. 

Finally, we caution that all our models adopted the 
on-the-spot approximation, in which recombination pho- 
tons are assumed to be reabsorbed in their zone of origin 
before escaping. In reality, reprocessed radiation isotrop- 
ically emitted within the ionized fingers may halt their 
growth by eroding them from within. Recombination 
radiation is believed to dominate the photon budget in 
the outer radii of centrally concentrated H II regions like 
those examined in our study (Ritzcrvcld 2005). Direct 
transport of diffuse ionizing radiation exhausts current 



ray tracing techniques, so the implementation of flux- 
limited diffusion for recombination photons in ZEUS-MP 
is currently underway in order to study their influence on 
instabilities in future models. Nevertheless, we expect 
that diffuse radiation will at most blunt the instabilities, 
not prevent their formation or destroy them because the 
fragmentation of the shell and breakout of the radiation 
through the cracks is a very robust feature of the fronts. 
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